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Abstract. Effective Poisson-Nernst-Planck (PNP) equations are derived for macroscopic ion 
transport in charged porous media. Homogenization analysis is performed for a two-component pe- 
riodic composite consisting of a dilute electrolyte continuum (described by standard PNP equations) 
and a continuous dielectric matrix, which is impermeable to the ions and carries a given surface 
charge. Three new features arise in the upscaled equations: (i) the effective ionic diffusivities and 
mobilities become tensors, related to the microstructure; (ii) the effective permittivity is also a tensor, 
depending on the electrolyte/matrix permittivity ratio and the ratio of the Debye screening length to 
mean pore size; and (iii) the surface charge per volume appears as a continuous "background charge 
density" . The coefficient tensors in the macroscopic PNP equations can be calculated from periodic 
reference cell problem, and several examples are considered. For an insulating solid matrix, all gra- 
dients are corrected by a single tortuosity tensor, and the Einstein relation holds at the macroscopic 
scale, which is not generally the case for a polarizable matrix. In the limit of thin double layers, 
Poisson's equation is replaced by macroscopic electroneutrality (balancing ionic and surface charges). 
The general form of the macroscopic PNP equations may also hold for concentrated solution theories, 
based on the local-density and mean-field approximations. These results have broad applicability to 
ion transport in porous electrodes, separators, membranes, ion-exchange resins, soils, porous rocks, 
and biological tissues. 
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1. Introduction. The theory of electrochemical transport in free solutions is 
well developed (57][65)[70] , but in many practical situations, ions move through porous 
microstructures with internal surface charge. Important examples in biology include 
nerve impulse propagation in the porous intracellular matrix of an axon 
ion transport through protein-based ion channels in cell membranes 



28 



89 



69 



selective 
and the 

electroporation of porous tissues for drug delivery and medical diagnostics 88 . In 
chemical engineering, the selective transport of ions and charged particles through 
membranes, gels and porous media is widely used for particle separations 
salination and ion exchange 



37 65 



71 



32 , de- 



energy 



characterization of porous rocks 
conversion in fuel cells [61] and energy storage in batteries [56] and electrochemical 
supercapacitors p5| . Analogous nanoscale transport phenomena are also beginning to 



be exploited in microfluidic devices 19 79 , which involve artificial porous structures 



with precisely controlled geometries and surface properties. 

In microscopic continuum models of electrolytes, the ionic fluxes are given by 
the Nernst-Planck equations describing diffusion and electromigration in the mean 
electric field, which is determined self-consistently from the mean ionic charge density 
via Poisson's equation. The resulting Poisson-Nernst-Planck (PNP) system has been 
studied extensively for the past century in the dilute solution approximation, not 
only for electrolytes [81112 29 , but also for semiconductors, where electrons and holes 
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behave like anions and cations, respectively [52]. The dilute-solution PNP equations 
can be derived rigorously from stochastic Langevin equations for the motion of point- 
like ions 



54 



Recently, a variety of modified PNP equations for concentrated solutions have 
been developed to describe strong interactions of finite-sized ions with charged surfaces 
at the nanoscale, as reviewed by J9|. Hard-sphere density functional theory 33 34 



and simpler mean-field models 43 62 have been used to modify the Nernst-Planck 

Poisson's equation has 
explicit 



equations for ionic fluxes to account for steric hindrance, 
also been modified to account for electrostatic correlations 



87 



11 21 35 72 



treatment of solvent dipoles 44 and solvation energy variations due to nonuniform 
permittivity 



All of these developments improve the microscopic description of 
ion transport close to charged surfaces, but our focus here is on the homogenization 
of such models over a charged microstructure to derive effective PNP equations valid 
at the macroscopic scale. 

There is a long history of heuristic models for macroscopic ion transport in charged 

A classical 



84 



membranes and porous media, dating back at least to the 1930s 
concept in membrane science, which we place on a rigorous footing below for gen- 
eral porous media, is the notion of a fixed "background charge" entering Poisson's 
equation, due to the volume-averaged surface charge of the porous medium [37| . In 
nanoporous membranes, the double layers are thick compared to the pore thickness, 
so that there are only small variations in diffuse ionic charge between the fixed surface 
or molecular charges. For most porous media, however, the double layers are assumed 
to be thin, leaving the pore spaces to be mostly filled with neutral solution, and Pois- 
son's equation is replaced by electroneutrality, without accounting for the background 
charge. In electrochemistry, this is a fundamental assumption of "porous electrode 
theory" (PET), introduced by Newman and Tobias |58 , which postulates electroneu- 
trality within the pores and effective Nernst-Planck equations of the same form as 
in the bulk solution, except for an empiricial tortuosity factor multiplying the ionic 
diffusivities. This approach has been applied extensively to batteries 26 45 56 57 . 
The nonlinear effects of double layer charging 12 have also recently been incorpo- 

The 



14 16 



rated into PET to model capacitive desalination and energy storage 
assumptions of PET have been tested against large-scale numerical solutions of the 
microscopic transport equations in certain cases of realistic microstructures 30 31 



but mathematical derivations are still be needed to predict the form of the macro- 
scopic equations and to provide a systematic method to calculate their coefficients. 
This is the goal of the present work. 

We are not aware of any prior attempts to rigorously homogenize the PNP equa- 
tions in charged porous microstructures, in spite of the many important applications 
listed above. Some formal derivations of PET have suggested such connections 14p5^ 



but without treating diffuse charge explicitly or relating macroscopic coefficients to 
the microstructure. For neutral species, the homogenization of linear diffusion over 
porous microstructures is well developed, and rigorous bounds are available for the 
effective macroscopic diffusivity tensor over all possible microstructures |85 . Homog- 
enization methods have also been applied to simplified models of ion transport, but 
not the full nonlinear PNP equations. Moyne and Murad [53] assume a Boltzmann 
equilibrium distribution of ions in a binary electrolyte at the pore scale and perform 
a homogenization analysis to derive effective equations for deformable porous media. 
Looker and Carnie |48 make the same approximation of microscopic Boltzmann equi- 
librium and derive symmetric Onsager relations for linear response, without stating 
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the general effective equations at the macroscopic scale. Allaire et al. [3] revisit the 
derivation of Looker and Carnie 48 using two-scale convergence methods developed 
by Nguetseng 59 and Allaire [l] (as we also use here) and prove the positive definite- 
ness of the Onsager tensor, but they make explicit use of an electroneutrality assump- 
tion. Although these studies simplify the treatment of diffuse charge in the pores, 
they do account for fluid flow coupled to ion transport. This leads to many complex 
electrokinetic phenomena in porous media and membranes |38j, w hose macroscopic 



mathematical description has been developed for decades 13 37 ,74 81 90 , albeit 
without the aid of rigorous homogenization methods that connect the microstructure 
to the macroscopic transport coefficients. For rigorous analytical results about the full 



Navier-Stokes-Nernst-Planck-Poisson system we refer the interested reader to 39 75 



and for according convergent finite element schemes to 66 



In this article, we derive porous-media PNP equations for charged microstruc- 



tures using two-scale homogenization methods. We extend the analysis of 77 78 
for the PNP-Stokes equations in uncharged dielectric porous media to account for 
the crucial, nonlinear influence of surface charge on the pore walls, although we ne- 
glect fluid flow. In the limit of thin double layers for isotropic media, our effective 



equations have the same classical form as those recently used by 27 36 50 68 , 82 



for nanochannels, where the potential is determined implicitly by macroscopic elec- 



troneutrality, including not only the ions, but also the surface charge. Helfferich 37 



credits Oel 60 with the formal justification of this approximation for charged porous 
media with thin double layers. Here, we derive more general PNP equations, valid 
for any double layer thickness, which preserve the form of Poisson's equation with a 
modified effective permittivity, where the electric field is produced by the total charge 
density. Our mathematical approach also provides a systematic means to calculate 
the tensorial coefficients in the macroscopic equations for a given microstructure. 

The article is organized as follows. We begin in Section [2] by recalling the PNP 
equations for homogeneous bulk solutions, and we summarize our results for porous- 
media PNP equations in section [3j We summarize the mathematical derivation using 
the two-scale convergence method in Section |4] In section [5j we briefly discuss the 
effective diffusivity and mobility tensors and investigate the validity of Einstein's 
relation between them. We illustrate our results for straight channels in Section |6] 
and irregular channels in Section [T] We discuss definitions of tortuosity in Section 
[8] and derive the general ambipolar diffusion equation for a binary electrolyte in a 
charged porous medium in section [9j following Mani and Bazant 49 . In Sections 



10 and 11 we take the limits of thin and thick double layers in the porous-media 
PNP equations, respectively. In Section |12[ we perform a simple microstructural 
optimization of the effective conductivity of a symmetric binary electrolyte for the 
case of parallel straight channels. In Section [TSj we conclude by discussing possible 
extensions and applications of our homogenized PNP equations. 

2. Poisson-Nernst-Planck equations for homogeneous media . 

2.1. Dilute solution theory . For simplicity, in this article we will perform 
most of our analysis for dilute, completely dissociated z : z electrolytes, focusing 
on their appropriate mathematical description in charged porous materials. We 
emphasize, however, that our results are by no means limited to symmetric bi- 
nary electrolytes and can be easily extended to multicomponent, asymmetric elec- 
trolytes. We adopt well studied mathematical framework for dilute binary elec- 
trolytes [9l[l2l[24[|42)|43l|62], which can also be used to describe general concentrated 
solutions (below) The concentrations of the ions c^{x^t) evolve according to mass 
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conservation laws 

dtc^ = -div (-c=^M±V/x±) , (2.1) 

where the classical Nernst-Planck fltixes (in parentheses) can be expressed in terms 
of the gradients of the ionic chemical potentials, given by dilute solution theory, 

H±= kTlnc^ + z±e(l). (2.2) 

For electrolytes (but not semiconductors) it is customary to neglect homogeneous gen- 
eration and recombination reactions, since charge transfer reactions tend to occur at 
interfaces. The variable (j) is the electrostatic potential, which describes the Coulomb 
interaction in a mean-field approximation. The coefficients D± are the (tracer) diffu- 
sivities of the two ionic species. The mobilities, M±, which give the drift velocity in 
response to an applied force, are then obtained by the Einstein relation M± = 
The total mean ionic charge density p controls the spatial variation of the potential 
(j) through Poisson's equation, 

-esA(j) = p:= ze{c+ -c~), (2.3) 

where e<j is the dielectric permittivity of the solution (roughly equal to that of the 
solvent), assumed to be a constant. The essence of the mean-field approximation is 
that each ion drifts in response to the electric field self-consistently produced by the 
mean charge density. For simplicity, we will mainly restrict ourselves to the case of a 
common diffusivity D = = £)_ and mobility M = M_|_ = M_ . 

In preparation for analyzing the full, nonlinear problem in a porous medium, we 
cast the equations in a dimensionless form using i as the reference length scale and 
to = £^/D as the reference time scale. We use the thermal voltage ^ as a scale for 
the electric potential. We introduce the reduced variables 

-4- c+ . c+ 7 zed) ^ X ~ t ^ /„ -s 

c+ = — , c- = — , <l>=^ x=- t=—, V = £V, (2.4) 
c c kl I In 

where c a reference concentration of ions, such as the nominal salt concentration 
of a quasi-neutral bulk electrolyte prior to its perfusion in the porous medium. The 
reference solution could be removed, or maintained in contact with the porous medium 
as a reservoir of ions at the reference concentration. We thus arrive at dimensionless 
Poisson-Nernst-Planck equations containing only the dimensionless parameter e = 

Ad 



diC+ = div + £+V(t)j , 

die- = div (Vc- - c-V<^) , (2-5) 
-e^A^ = c+ - c~ , 

where Ad is the Debye screening length of the reference solution, defined by Xd ■= 

1 /2 

^'-'^ ^ symmetric binary electrolyte. 
In the subsequent sections, we also work with the dimensionless charge density p 
and total ion concentration c, which are defined by 



. c+ + c" 
c= — ^= — > P = 



2c ' " 2c 



(2.6) 
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respectively 8]|T2 . In terms of these dimensionless variables (2.6), the PNP equations 



take an elegant symmetric form, 

dfc ^ div (Vc + pV(f^ , 

dfp = div (V/5 + 5V^) , (2-7) 



-e^A4> ^ p. 

This formulation is not only mathematically convenient, but also provides physical 
insight by distinguishing between dynamical phenomena involving diffuse charge, such 
as capacitive charge storage and surface conduction, and those involving neutral salt 
transport, such as bulk diffusion, surface diffusion and capacitive desalination |12[|14l 



23 24 62 . This formulation is also useful in the analysis of Frumkin double-layer 
effects on Faradaic electron-transfer reaction kinetics (8| [T6{[T7) . 

In our analysis below, we shall use dimensionless equations and drop the tilde 
accents for ease of notation. 

2.2. Concentrated solution theories . In a concentrated solution, various 
physical effects contribute an excess chemical potential of an ion, /x^^, beyond that of 



a dilute solution (2.2 1, 



p, = kTln c' + z,e0 + p'^'= = kTln{j,c') + z^ecf) (2.8) 

where we also define the molar activity coefficient 7^ = exp(/i^^/fcr), and the diffu- 
sivity and mobility may also contain nonlinear modifications [o]. The excess chemical 
potential is generally a non-local functional of the ion concentrations and potential, 
which can be treated using statistical density functional theory |34| or the weighted- 
density approximation [4^, leading to integro-differential equations. In the local- 
density approximation (LDA), the excess chemical potential depends only pointwise 
on the ion concentrations, and the model leads to partial differential equations. The 
LDA becomes exact in a uniform bulk solution, but it is also widely used in the pres- 
ence of strong concentration gradients, although not always maintaining accuracy. 



Classical concentrated solution theory 57 assumes a neutral bulk solution, and the 



concentration-dependence of the activity coefficient and diffusivity are left for empir- 
ical fitting to experimental data. Alternatively, concentrated solution models can be 
derived from various microscopic physical assumptions, as in charged hard-sphere or 
lattice-gas models |;9]. 

3. Homogenized Poisson-Nernst-Planck equations for porous media. 

3.1. Microscopic model. Our derivation begins with a periodic representation 
of the microstructure and uses a formal two-scale convergence method to derive effec- 
tive equations valid at scales much larger than the geometrical period, as sketched in 



Figure 3.1 Schmuck 77 recently apphed this approach to derive effective PNP equa- 
tions for ion transport with convection in uncharged porous media. The extension 
to fluid flow for the full nonlinear PNP equations in a charged microstructure is still 
lacking a consistent derivation, as pointed out in 76 . Indeed, the classical theory of 
electrokinetic phenomena in porous media (e.g. streaming potential, electro-osmotic 
flow, diffusio-osmosis) assumes linear response to small perturbations [38| , but it has 
only recently been understood that many nonlinear electrokinetic phenomena, such as 



second-kind electro-osmotic flows 46 83 , induced-charge electro-osmotic flows 10 
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deionization shocks 49 50 , and surface-driven over-limiting current 27 , can arise 
due to the couphng between ion transport and surface charge in microstructures. 
Hence, we neglect flow, but we allow for a fixed surface charge density on the solid or 
molecular matrix. In the following sections, we emphasize physical implications and 
provide some illustrative examples. 

Suppose that the heterogeneities, i.e. the size of the pores defined by a reference 



cell Y, as for example in Figure 6.1 left (in 2D), are very small with respect to the 
size of the porous medium denoted byfJcM'', l<(i<3, and that they are evenly 
distributed. This is a realistic assumption for a large class of applications. From a 
mathematical point of view, one can model this distribution by supposing that it is 
a periodic one. In fact, to reduce the approximation error, one can first compute the 
averaged pore structure of a suitable periodic replacement by a so-called Y- Algorithm, 
The periodicity can be represented by a small parameter r and hence the 



77 



see 

porous subset il^ of is then denoted by Qf. and correspondingly the solid subset 
ri'' by ri* := f2 \ f2^. For the moment, we assume that the solid-electrolyte interface 
Ir :— dflP nd^l^ is smooth. Further, we suppose that we are given a reference period 
Y, in which the reference heterogeneities are defined, see Figure 6.1 left (in 2D) 



respectively Figure 6.1 right (in 3D). The pores in Q are periodic of period rY and 
their size is of order r. The concentration densities have to be solved only in the 
domain of the liquid phase which is the perforated domain fl^ = i}\nf.. Then, for 
the scaling parameter r > 0, the positive ion density , negative ion density , and 
the electrostatic potential (j)r are governed by the equations 



(micro PNP) 



dtc+ = div (Vc+ + c+W(pr) in , 

inQP, (3.1) 

in n:=npuni, 



dtc,, = div (Vc^ — V0r) 
—div {e{x/r)\/(f)r) — — 



where e{x) 



the dimensionless Debye length, and 



the dimensionless electric permittivity of the solid material Eg scaled to that 



of the pore phase e^. The system (3.1 1 is completed by the boundary conditions 



\I^(j)r = onir := dnP n d^l 



VnCr - C^Vn(t>r =0 On , 

-£(x/r)V„0r — r<Js{x/r) on 1^ , 



(3.2) 



where <Ts{x/r) is a periodic surface charge density. Since we describe the solid- 
electrolyte composite by macroscopic respectively continuous variables, we do not 
model the surface charge by a jump boundary condition. In Figure [3?l] we illustrate 



the periodic formulation (3.lH(3.2), which allows to introduce a periodic reference 



cell Y := YP U Y'^ capturing the pore geometry. 

3.2. Macroscopic approximation . A formal application of the two-scale con- 
vergence allows to derive the following upscaled equations describing electrolytes in 
porous media for fixed surface charge. The basic idea is an asymptotic approximation 
of the form 



1 



0{r'))4{t,x) 



(3.3) 



Mt, x)^{l + r6'^{x/r)d^, + 0{r^)) M^, x) , 
where Einstein's summation convention is applied over repeated indices and the mi 



croscopic variable y = x/r e F is used. We point out that the series in (3.3) are 
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Macro scale: (Xi, X2) 



Micro scale: (y^, yj 




Figure 3.1. Left, macro scale: Domain Q := U Qf. with the solid-liquid interface Ir := 
dO.^ n aO^. Right, micro scale: Reference cell Y := YP VJ Y^ := [0,^i] X [0,^2] with solid-liquid 
interface S := Y^ D Y" = Ir f\Y . The micro scale is obtained by rescaling the macro scale with 
the parameter r which measures the characteristic size of the heterogeneities. The upscaling or 
averaging process then consists in passing to the limit r — > 0, i.e. the electrolyte and the solid phase 
are homogeneously mixed under keeping the corresponding volume fractions constant. 



not necessarily convergent. By Y we denote the reference cell which captures the 
geometric and interfacial information of the microstructure. 



It is important to stress that we use dimensionless length and time variables (2.4 1 



scaled to the characteristic length and time scales for diffusion across a reference cell 
at the scale of individual pores. This is convenient for homogenization and allows us 
to vary the Debye length relative to the mean pore size (below). In applications of 
the macroscopic PNP equations obtained by homogenization, however, it is better to 
rescale length and time according to the scales of bulk diffusion across ^ 1 refer- 
ence cells over a length L — (./r uv terms of the following macroscopic dimensionless 
variables: 



L 



LV 



^V, t 



tD 
L2 



(3.4) 



This rescaling will be useful below when we consider the limit of thick double layers 
compared to the pore size (e ^ 1), which could be thin (e ^ 1) or thick (e ^ 1) at 
the macroscopic geometrical scale. 

Without further details, we directly state here the main result which extends the 



derivations in 77[ 78 to the case of charged porous materials. We further generalize 

The article 



77 under local thermodynamic equilibrium as achieved in 76 78 



78 



additionally qualitatively characterizes the effective porous media approximation of 
the PNP equations by quantitative error estimates. 

The following definition guarantees the well-posedness of the effective porous me- 
dia approximation subsequently presented. It can also be justified on both physical 
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and mathematical grounds by the assumed separation of the macroscopic and micro- 
scopic length scales, analogous to the quasi-equilibrium approximation for thin double 
layers. 

Definition 1. (Local thermodynamic equilibrium in dilute solutions) We say 
that for dilute solutions the reference cells Y are in local thermodynamic equilibrium 
if and only if, 

fiQ := log c^tjx) + Zi(j>{t,x) = const. inx/s = y€Y, (3-5) 

where (Iq denotes a constant value of the chemical potential of positive (i = +) and 
negative (i = —) ion densities and z± := ±1. The locally constant potential /Iq can 
assume different values in different reference cells. 

Remark 1. We emphasize that the equilibrium condition (3.5) specifically ac- 
counts for dilute solutions. For the more general concentrated solution theory, a cor- 
responding equilibrium is defined in \2.^ . 

Under local thermodynamic equilibrium, the effective positive ion density c"*", 
effective concentration of negative ions c^, and the effective electric potential (f) are 
determined by 

' pdtc+ = div (t)Vc+ + c+MV^] , in 17 , 

(macro PNP) <| p^^c- = div mVc" - c-MV(/)j , in f7 , (3.6) 

-div (£(e, q;)V(?!)) — p (c+ — c~) + Ps , m Q , 

where Ps •= jyy IdYindY2 '^^ ^^iv) total surface charge per volume (which plays 

the role of the "background charge" in classical membrane and semiconductor mod- 
els), a := 1^ is the permittivity ratio, e :— \J 2 e / ^ dimensionless De- 
bye length (measuring the thickness of the double layers), and the different tensors 
D := {<^ki}i<k.i<di M {mfc/}i<fc ;<^, and £(e, a) := {e«(e, a)}i<fe,,<<i are defined 

by 

dfc/ = 1^ {^ki- Skjdy^C' {y)} dy for i = f , 2 , 
nife, = 1^ {^ki- Skjdy^^^^' iv)} dy for i = 1, 2 , 

eki{e,a) ^ ^J^^{Ski - SkA^e^'iy)} dy + ^^{Sm - Sk,dy^e''{v)} dy . 

(3.7) 

The corrector functions ^""^ solve for i — 1,2,3 and 1 < r < d the reference cell 
problems 

-Ayie^ -yr) = -Ay{e'''- -yr) 

(_V^(e'^ - yr) + ^yie^^ - yr)) ■ nj ^ ou I := dY' , 

^iv is ys_periodic and My^ (C"'') = , (3.8) 
-diyy{e{y)Vyie^'- ^yr)) =0 in F , 

^33^ is F-periodic and M {^^^'' ) — ■ 



As a consequence of (3.8 1, it holds — ^22,. _ rpj^g quantity dY'^ denotes the 

boundary of the solid phase in the reference cell Y = Y'^ U Y^ and is part of the 
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interface / on the domain il, see for example Figure [6?T] left. The subset denotes 



the solid phase (black in Figure 3.11 and the electrolyte phase (blue in Figure 
3.1). Furthermore, the effective tensors (3.7li-(3.7)2 contain relevant macroscopic 



corrections defined by the microsc opic pore structure. We remark that the corrector 
£fc/(e,a) for the Poisson equation (3.6)3 is equal to D^i defined in (3.7)i, if and only 
if we assume that the electric field only exists in the electrolyte phase. Such an 
assumption is equivalent to the limit a — ^ 0, see Section (3.6 1. The correction tensor 
Dfc; can only be motivated analytically in the case of straight channels, see Section |6] 
and [sj. Already for perturbed straight channels, a numerical approach is required to 
solve the reference cell problem (3.8). We refer to Section[7]and for further details. 

Finally, we motivate that the assumption of local thermodynamic equilibrium 
(3.5) gu aran tees solvability of the reference cell problem ( |3.8[ ). As demonstrated 
in |76l, (3.8)i is obtained as the C'(r~^)-order problem for an asymptotic expansion 

: 1,2, we first obtain 



That means, for arbitrary smooth test functions 'J* 



k,l=l ^ 

+ {{c'Skidy,e^'-,dy,^')^^ - {c'5kr,dy,^^)^^) zA^cj^^ =0 for * = 1,2, 

(3.9) 

where = c"*", = c~, zi — +1, and zi — —1. The notation (u,v)d ■= Jj^uvdx 
denotes the scalar product in L^(£')-sense with D = Y''. In classical homogenization 
theory, one can immediately cancel out macroscopic quantities like d^^d ■ This is here 
not possible due to the additional gradient of the electric potential d^^ip- Since we 
do not exactly know when dx^c^{x,t) becomes zero, we cannot immediately apply 
Lax-Milgram's theorem or the Fredholm alternative in order to guarantee solvability 



of (3.9 1. But under the local equilibrium assumption (3.5), we are able to remove the 



macroscopic quantities d^^c^ and Zidx^4> such that we end up with equation (3.^ 
Solvability of (3.8)i then immediately follows via Lax-Milgram's theorem. 



3.3. Einstein relations. The upscaled PNP equations indicate that Einstein's 
relation between diffusion D and mobility M coefficient, i.e., M ~ does not hold 

with respect to the porous media correction tensors D and M. At first, this seems 



to be physically inconsistent, since we seem to loose the gradient flow structure (2.1 1 
and a Boltzmann distribution for ion densities in equilibrium. 

However, we point out that the tensors D and M are corrections to Vc^ and V0, 
respectively, and not to diffusion and mobility coefficients. This important fact then 
motivates to define the mean field approximations 



Vc± := DVc± 
V0 := MV(/). 



(3.10) 



It is immediately clear that the Einstein relation holds for the quantities (3.10 ), i.e 



D±\/c^ + kTz±M±V(j) — fcTM±V/i_|_ in dimensional variables. Moreover, via (3.10) 
we can define a mean field gradient of the chemical potential by 
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such that gradient flow (2.1) becomes 



"div I 



(3.12) 



where M — M± — We further remark that the mean field approximation (3.111 



is only defined for the gradient of the chemical potential. This fact is a direct con- 



sequence of the Ansatz (3.3). As a consequence, the chemical potential /i± remains 



unchanged and hence Boltzmann's distribution for the ion densities still holds in ther- 
modynamic equilibrium. 

3.4. Material tensor and Onsager relations. In [77 , it is shown that the 
definitions ( 3.7 )i-( 3.7)3 allow to introduce a so-called effective material tensor 













which allows to write (3.6) for the field vector Q 
side I(Q) := [0, 0, -Q^' by 



:= [c+,c- 



5tQ + div (s(Q)Vq) =I(Q), 



where dt is the operator 



dt 










dt 




(3.13) 

]' and the right-hand 
(3.14) 

(3.15) 



and also V and div are correspondingly defined. 

Remark 2. Let us point out the special case where the electric field only exists in 
the electrolyte phase. In su ch situations, the tensor iki{e, a) in (3.13) can be replaced 
by e^Dfc/, see Section (5, 



The material tensor ( |3.13 1 provides a "nonlinear equivalent" to the Onsager re- 
lations for situations apart from thermodynamic equilibrium. It is completely deter- 
mined by the elementary cell Y as for example given in Figure 6.1 left (in 2D) and 



right (in 3D), respectively the equations (3.7)-(3.8). Such a reference cell represents 
the characteristic geometry of the porous material under consideration. The mate- 
rial tensor S relates the gradients of concentrations and electric potential to their 
corresponding fiuxes J := [Jc+ , Jc- , J^]', i-e. 



J = S(Q)VQ. 



(3.16) 



The classical Onsager reciprocal relations only hold true for a linearized material 
tensor (3.13). In fact, Allaire et al. [3] rigorously achieve the Onsager relations, first 
verified in Looker and Carnie |48| , for a linearized and time- independent version of 
ionic transport in contrast to the full nonlinear problem considered here. 

Next to the physical information gained by the transport parameters in the tensor 



(3.13), the effective equations (3.6) also prevent high-dimensional and hence compu- 



tationally expensive problems. Such computationally demanding situations gener- 
ally result by solving (2.5) over the real microstructure. In fact, the mesh for such 



computations is required to be smaller than the characteristic length scale r of the 
heterogeneities. 
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Let us mention that it is not possible to derive (3.6) by volume averaging or the 



representative volume mehtod (RVM), since the system (2.51 is nonlinear. Moreover, 



it is not clear what is the right size of the test volumes with respect to which the 
gradients and their corresponding fluxes are averaged. Also a possible dependence 
on a source term and boundary conditions, see (3.2)3, cannot be captured by such 



approaches. These reasons strongly motivate the effective equations (3.6) derived here 
by homogenization techniques. 



3.5. Reformulation of (3.6) for salt and charge variables. Finally, we 



reformulate (3.6) for the physical quantities charge p and salt c, i.e. 



pdtp = div {t)Vp + cMV0^ , in 17 , 

pdtc = div {hVc + pMV0^ , in f2 , 
—div {e{t, a)V(f>) = pp + ps , in fl . 



(3.17) 



where we call k :— pM the macroscopic conductivity. D and M are microscopic 
correction tensors which contain the information about the pore geometry. 

There is a growing body of recent theoretical and experimental work on ion trans- 
port in charged nanopores. If the Debye length is smaller than the pore radius, then 
continuum models can agree very well with experimental data 68] . The mathematical 



description of ion transport in various recent models 27 49 50 



67 82 has a similar 



form as our porous-media PNP system (3.171, but without making any connection 
between the microstructure and the macroscopic coefficients in the equations. In par- 
ticular, the hydrodynamic factors 



and Kj r used in 



82 



are not systematically 

derived and not defined by the pore geometry as D and M in ti.7\ and (3.8). More- 



over, the cross-sectional area averaging developed by |67_ and [50] was only applied 
to straight, symmetric channels. 

3.6. Insulating porous matrix. As indicated in Remark [2] the porous media 
correction e(e, a) for the Poisson equation can be reduced to D in the case of an 
insulating porous material. As a consequence, also M reduces to D and hence the 



material tensor (3.13) simplifies to 



Dfci Q'^'Dki 
Dm -QiDfe, 
Ski 



(3.18) 



Let us introduce the following coordinate tranformation 



D 



-1/2, 



(3.19) 



where components of x admitting "cx)" are subsequently to be treated as parameters. 
This coordinate transformation can also be interpreted from a more experimental 
view point, since there is a close relation to tortuosity and diffusivity, see Section [8] 
for details. With (3.19) the gradient Va, and the divergence operator div^: change 



with respect to the new coordinates to 



= D^^^Vs , and div^ (V^^)' ^ divjD"^/^ . 



(3.20) 
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Via (3.20), the tensor (3.18) can be written in this new coordinates x by 



I 



(3.21) 



where Q(t,I)^/^x) — Q{t,x). Hence, the material tensor (3.18) takes the same form 
in the new coordinates x as the classical PNP eqations for homogeneous media in 



the case of an insulating porous matrix. Moreover, the porous media equation (3.141 
reads in the new coordinates as 



dtQ + divi (M(Q)ViQ) = I(Q) . 



(3.22) 



Let us remark that the transformation (3.19 ) accounts for a finite separation of scales 



and can be generalized to the case of a continuum of scales by the idea of metric-based 
upscaling introduce in . 

3.7. Concentrated solution theory. We briefiy consider arbitrary chemical 



potentials (2.8 1 satisfying the local density approximation, i.e. depending only on 



local concentrations and the potential. In that case, the gradient flow (2.1 1 can be 
accordingly obtained for species and electric potential 4>. To this end, we first 



compute the gradient of the chemical potential (2.2 1, i.e 



V/x(c' 



(5c* ~^ 66 ^ 



(3.23) 



If the first variations |^ and ^ are linear, then the form of the upscaled equations 

the multiple-scale method 



can be verified as in Section |4j For nonlinear 



1^ and 



(3.3) allows us to at least formally obtain the structure of the homogenized equations. 



In light of the explanations after (3.23), we expect that the general porous media 



formulation for LDA models takes the same basic form 



pdtCi := -div M^V ^{cj , (f))) 



(3.24) 



where V/i is defined as in (3.11) and we allow for a multicomponent electrolyte, as 



well as (formally) for nonlinearities in the mobilities. 

4. Formal derivation of the upscaled equations by the two-scale con- 
vergence method. In Figure [4?T] we recall the basic idea behind the homogenization 
method. For the upscaling process, we first write the system (3.1) in the distribu- 



tional sense. Therefore, we multiply equations (3.1 )i and (3.1)2 by a smooth function 



if G C°°(r2^) and the Poisson equation (3.1 )3 with ip G C°°{Q). After integration over 
and f2, we end up with the formulation 

(9tC+,V3)j^p = (Vc+ -f CrV0r,Vv3)j^p + J {V pr + CrV (t)r} do{x) , 

{dtc~,ip)^p = (Vc7 +c+V0r, V(/3)j2P + J {Vc~ + c:^y (j)r} ipn do{x) , 
- {e{x/r)V(j)r,VLp)^-^ = (c,+ - c~,(p)^- J e{x/r)V4>ripndo{x) . 

(4.1) 
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Strongly heterogeneous 
material, i.e. < r ^ 1 



□nnnnn 
□□□□□□ 

□□□□□□ 
□□□□□□ 



Homogenous approximation 



(r ^ 0) 



Figure 4.1. Left: A composite material whose characteristic heterogeneity has the lenth r. 
Middle: Passing to the limit r — > under constant volume fraction between circle and square. 
Right: The limit problem is obtained with the theory of homogenization. 



With the boundary conditions (3.2 1, the system ( |4.l[ ) reduces to 



(4.2) 



- (e(x/r)V0r, Vv3)fj = (c+ - + r J <Jsipdo[x), 

where the boundary conditions on need to be scaled by r as motivated in [2| |55| . 
Now, we take the limit r — )> and get the upscaled or homogenized description of 
(4.2). We recall that the concentrations need to be extended to the whole domain 
Q. in (4.2) before one passes to the limit r 0, as in 77 for the case without surface 
charge. Hence, it remains to account for the presence of the periodic, surface charge 
crs{x/r) on the boundary 1^ :— dil^ O dilf.. Therefore, we apply the test function 
iprix) := <y3(a;) + ripi(x^ ^) to (4.2)3. Subsequently, we only consider the boundary 
integral, i.e.. 



r / as{x/r)ipr do{x) 



1 

W\ 



(p{x) 



dYindY2 



fTs(y) dxdo{y) 



(4.3) 



Using (4.3) and the limit in the other terms of (4.2)3, achieved in [TT] (a nd in a more 
general context in ^76]), finally leads to upscaled Poisson equation (|3.6[)3. 



5. Effective diffusivity and mobility. In order to get the effective diffusivity 
tensor D and mobility tensor M, we need to go back to the dimensional quantities C±, 
which represent concentrations of charged species in the considered electrolyte. Since 
we restrict ourselves to a binary symmetric electrolytes, we only distinguish between 
positive charged species C+ and negative ones C_. In these variables, the upscaled 
equations look like 



pdtC± 
-div (eo(ep,a)V<I>) 

where we used the definitions 



: div (D±VC± ±zeC±I 
: pze (C+ - C_) + ps , 



zV$) 



J± :=i?±D, 
4I± A//±M . 



(5.1) 



(5.2) 



14 



M. Schmuck and M. Z. Bazant 



We recall that the mobility is defined via Einstein's relation by AI± — Hence, the 
relations (5.2) define the dimensional effective diffusivity tensor D± and the dimen- 



sional effective mobility tensor M±. Furthermore, the equations (5.2) further indicate 
that (Nernst-)Einstein's relation is still consistent with the upscaled formulation (5.1 ). 
However, since the tensors D and M are corrections to the gradients and not to D± 
and M±, Einstein's relations do not hold for the effective quantities D± and M± 
see Section [3| ( |3.3[ ). We also emphasize that the applicability of the Nernst-Einstein 
equation must not be guaranteed in every situation, see for example the discussion 
in 57 Section 11.7 and Section 12.5]. There, it is pointed out that the validity of 
Einstein's relation is valid at infinite dilution, although its failure is related to the ap- 
proximate nature of the flux of the Poisson-Nernst-Planck equations. But primarily, 
the validity of the Nernst-Einstein relation rests on the fact that the driving force for 
both migration and diffusion is the gradient of the electrochemical potential, and the 
decomposition of this into a concentration term and an electrostatic-potential term is 
without basic physical significance. 

6. Solving the reference cell problem for straight channels. For the sub- 
sequently studied examples, wc additionally suppose that the electric potential </> only 
exists in the electrolyte phase as the salt and charge concentration. This assum ption 
implies that the correction tensor e(e, a) becomes eDki as motivated in Section|3|(|3.6[). 



The two-dimensional case: For simplicity, let us assume the molecular diffusion 



tensor D as isotropic, Di 



Hence, via Einstein's relation, i.e. Mi. 



^ the 



mobility tensor M is also isotropic. We consider a reference cell defined by Figure 6.1 
left (in 2D). The porous media correction with respect to the diffusion can be written 
in the two-dimensional case as follows 



D 



Dii 




Oviously, in the case considered we have, as in '5', Dn = pD and D22 — 0. Hence, the 




D22 



(6.1) 



porous media Poisson-Nernst-Planck equations (3.6) take the following elegant form 

t^c + d.,,{pdx^c^) inf], (6.2) 
'pedl^<p = pp + Ps inn. 



dtc 



The three-dimensional case: The porous media correction tensor for the diffusion 
is obtained by the same arguments as in the two-dimensional case, i.e. 



(6.3) 





" Du 








D = 





D22 













D33 . 



where Du — D33 — p and D22 = 0. A strightforward extension of the straight channel 



to dimension three is depicted in Figure 6.1 right. As in (6.2), we obtain the following 



porous media equations for straight channels in the three-dimensional case, i.e. 

dtp = {dl^ + dljp + (9., {cd^,4>) + (ca,3 0)) in n , 

{dl C) c + {d., (p9.,0) + ipd.A)) in ^ , (6.4) 



dtc 



-pe {dl^+dljc^ = pp + ps 



in n , 
in . 
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Figure 6.1. Example of straight channels: Left: Two-dimensional case, (pore phase is 
red) Right; Three-dimensional case. 




Figure 7.1. Perturbed straight channels in 3D, see |6,: Left: Reference cell geometry. 
Right: Cross-section of the period. 



7. Solving the reference cell problem for perturbed straight channels 
(3D). In difference to straight channels, the case of perturbed straight channels re- 
quires the numerical calculation of the components da for i = 1, 3 of the effective 
diffusion tensor D. However, the component is as in (6.3). For details about the 
numerics, we refer to [6] and state immediately the computational results here, i.e. 



D 



0.3833 

1 



(7.1) 



Hence, with (7.1) the porous media approximation (6.4) immediately changes in case 



of perturbed channels as depicted in Figure |7.1| to the following system of equations 
dtp - (0.383391^ + dlj p + (d,, (0.3833c9,,(/.) + 8^, {cd,,<j))) in n , 



dtc = (0.38339^^ + dlj c + id,, (O.3833pa,,0) + d,, [pd^M in 17 , 

-pe (0.38339^^ + 5^ J = pp + in 17 . 



(7.2) 
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Let US now compare system (7.2) for perturbed straight channels with the straight 
channel example (6.4). One immediately recognizes that the porosity parameter p 
doesn't affect the ion density equations. Further, it is apparent that longer fluid 
paths modify these ion concentration equations accordingly. This fact is connected to 
the study of the tortuosity, which is the topic of the next section. 

8. Discussion of tortuosity and efTective diffusivity. In the following, we 
motivate that homogenization allows to validate current tortuosity relations and to 
give directions towards refinements of such relations. The explicit examples from the 
Sections [6] and [7] allow to systematically understand the influence of the geometric 
structure to the tortuosity. Sometimes, people introduce a so-called diffusibility Q 
which relates the molecular diffusion constant D{ to the effective diffusion constant 
Dp of a porous medium, i.e. 

Dp = QDf . (8.1) 

The expressions for Q available in literature can be divided into three classes, see 
Brakel et al. [86]: (1) empirical correlations, (2) semi-empirical equations based on a 
pore model, and (3) theoretical expressions. 

(1) The empirical correlations express Q as a function of the porosity p, i.e. 
Q — f{p). (2) If Q is defined by the special class of functions f{p) = jp^. The term 
p^ is generally said to account for the influence of the smaller cross sectional surface 
available for diffusion. (3) The theoretical expressions for Q have been derived for 
dispersed solids in the form of spheres. 

But, let us first give a historical overview. In any porous system, the presence of 
solid particles/material causes the diffusion paths of species to deviate from straight 
lines. Consequently, the diffusion coefficients of species must be corrected. One tries 
to capture this deviation from straight lines in a porous medium by a term called 
tortuosity, whose suitable definition is still an actual research topic. This is also 
motivated by the following considerations. 

By theory and dimensional reasoning, Petersen [64 suggests that the diffusion 
coefficient is scaled with tortuosity as follows 

Dp^^, (8.2) 

which implies for the diffusibility Q = 1/r^. A similar relationship is introduced by 



Aris j4j and Satterfield [73], i.e 



Dp = ^Di , (8.3) 

T 

and hence Q = p/r. The simplest and most intuitive method to estimate the tortu- 
osity T (in the one-dimensional case) is the ratio of the length of the curve to the 
shortest distance of its end points Lab, i-e. 

r:=^. (8.4) 

^ab 



In Brakel et al. 86 , they consider a slight generalization of (8.2) by introducing a 
constrictivity parameter d := (^^^^ , which accounts for the fact that the cross 
section of a segment varies over its length. Hence, (8.2) changes to 



(8.5) 
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sucht that Q — 

Further, Brakel et al. 
Q — f{p) does not exist. Moreover, they emphasize that the pragmatic value of the 
available Q — p relations is not very great. 



86 argue that for porous materials a function of the kind 



80 



give a 



Recently, also Shen and Chen 
critical review of the current impact of tortuosity on diffusion. Therefore, we motivate 
our discussion and study of Q in this section by suggesting a theoretically obtained Q 
with the help of homogenization theory. The diffusivity Q could turn out as a relevant 
parameter to compare empirical measurements with theoretically obtained effective 
quantities. 

Up to this end, we first extend the above relations to tensorial versions, i.e. we 
denote by Dp the effective diffusion tensor in a porous environment and by Df := 
{DSij} ■ the molecular diffusion tensor in free space. 5ij denotes the Kronecker delta 



function. First, we extend (8.2) to 




1/2 



(8.6) 



where Dp :— DD and the diffusion corrector D is obtained by homogenization in 



(3.7)i. The quotient in (8.6 1 is understood as a division of each component if and only 



if both, the corresponding nominator Df and the denominator Dp are different from 
zero. Instead of this restriction, we could also allow for a tortuosity f with components 
admitting "oo" . We point out that the tensorial relation 
diffusivity, i.e. Q = l/f^. 

Another very interesting interpretation of 



.6 ) also implies a tensorial 



sulating porous matrix, see Section 3.6 The tortuosity f in 
coordinate transformation (3.19), i.e. x = tx. 
In view of 



6|) is possible in the case of an in- 
6 ) corresponds to the 



.3) and (8.5), we motivate further the extensions of (8.3) to 



pDj 

Dn 



■7) 



with corresponding Q = p/t and 



pdjjDf 

Dr. 



1/2 



with Q = pd/f^ 



Let us apply definition (8.6 1 to the examples from Sections [6] and [T] In the case 



of straight channels, see Figure|6.1|on the right-hand side, the definition (|8.6|) implies 
the following tortuosity tensor 







i/Vp 



(8.9) 



We point out that the porosity p with respect to straight channels corresponds to 
the channel height on the unit reference cell. Let us compare (8.9) with the intuitive 
definition (8.4). If we apply definition (8.4) in a straightforward manner, then t = 1. 
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However, it is not clear for straight channels, which path L-y is reasonable, 
check for example the average 



1 ^ 

i=l 



Let us 



(8.10) 



With 



.10), we get r 

Lab — 1 



L-y where Lj — (^\/l + + 1 + (1 + p)^ for the path 



^7 

lengths 7i Lab — 1, 72 := Lab+P = 1+p where the porosity p is the channel height, 
and 73 := \/l +p'^ the diameter. But due to Boudreau 18 Section 2], the tortuosity 



must approach unity fo rp -> 1. This is violated b y defi nition ( 8.4 ) together with ( 8.10 1 

But the tortuosity (8.9) defined via the homogenization 



1 



V2 



with 



.10 



process perfectly satisfies this condition, see [Isj Section 2] . Accordingly, in the case 
of perturbed straight channels as considered in Section [t] the tortuosity tensor 
becomes 



.61 



1/V0.3833J5 


1/^ 



(8.11) 



One immediately recognizes that fn in (8.11) is > 1 in the limit p — ;> 1. Hence 



Boudreau 18 Section 2] doesn't hold. This contradictions require caution by using 
definitions (8.2) and (18. 61). Next, we examine the definition (8.7) which becomes for 



the case of straight channels 



1 

1 



A comparison of 

However, in the case of perturbed straight channels 



12 ) with ( 8.4 ) shows perfect agreement, i.e. (rigj^ 



11 



.12) 



we depend on the numerical 
accuracy. Since the mesh in |6 is not very fine, we cannot necessarily expect equality. 
In fact, we obtain (r^ jg"]^ ) — 2.6 and r jgTj^ = ^^^^ = 1-9. However, this discrepancy 
also motivates the critical statements of 80 , 86 about the pragmatic value of tortuosity 



as mentioned above. We leave the investigation of the definition ( 8.8 ) to the interested 
reader, since the definition of the constrictivity parameter in 86 is a delicate point 
and again a new source for modeling errors. 

As a conclusion of this discussion, we motivate that homogenization theory al- 
lows to derive effective equations which do not require a questionable tortuosity or 
diffusivity parameter. Moreover, the correction tensors obtained by the two-scale con- 
vergence method provide a tool to check available tortuosity or diffusivity definitions 
and might provide directions how to improve their consistency. 

9. Ambipolar diffusion equation for a binary electrolyte. Motivated by 
the considerations in Mani and Bazant [49] by volume-averaging, we study here the 
same problem with the help of homogenization theory. The advantage relies on the 
fact that we are able to accurately treat the nonlinear terms. Up to now, there exists 
no general rule how to upscale the nonlinear terms by volume-averaging approaches. 



Hence, we extend the porous media approximation (3.6) to a dilute, asymmetric 
binary electrolyte {i ~ +, — ) with arbitrary ionic charges, q± = ±z±e in this section. 



Ion transport equations in porous media 



19 



For simplicity, we assume constant diffusivities D± in the microstructure and denote 
the corresponding upscaled diffusivities and mobilities by D± = D±t) and M± = 
M±D. Without loss of generality, we consider a negative surface charge, i.e., ps < 0. 
Moreover, we work in the context of an insulating porous matri x such that the porous 
media correction tensors satisfy D = M = e, see Section [sj |3.6| . We simplify now the 
Poisson-Nernst-Planck system by applying the usual conventions 



= pe(z+C+ - z_C_) + 
pc — p (z+C+ + z_C_) + 



(9.1) 



where the first relation expresses quasi-neutrality for the case of surface charge. This 



assumption naturally arises here in view of the derived effective equations (3.17) for 
fixed surface charge. However, in [51 67 82 such a neutrality condition has been 



suggested by purely physical reasoning. Furthermore, we will not make use of the 
Nernst-Einstein equation between the diffusivity tensors D± — D±D and mobility 
tensors M± — M±D. Hence, the ambipolar diffusion equation derived under such 
assumptions takes the form 



pdtc = i:»div (dVc) 



- -div (^PsDV^ 



D.z 



kTez^M^ 



-div 



where we used relations 

z+M+D_ + z-M^D+ 



D := 



Z+M+ + z_M_ 



and 



2z+z_M+M_kT 
z+D_M+ + z_D+M_ 



(9.2) 



(9.3) 



We remember that D is defined as (3.7)i. 



The correction tensors D for straight channels (in Section |6| and for perturbed 
straight channels (in Section [t]) allow to accordingly rewrite the ambipolar diffusion 
equation (9.2), which describes a porous material for a surface charge density Gs- In 



view of the volume-averaged straight channels studied in 49 , we only consider in the 
following the example from Section |6] With (6.1), the equation (9.2) immediately 
takes the form 



dtc : 



D. 



kTez^M^ 



(9.4) 



An interesting remark regarding (9.4 1 is that the porosity parameter p cancels out. 



10. Thin double layers at the pore scale . Recently, Mani et al. J50\ devel- 
oped a " simple model" describing the propagation of salt gradients in a microchannel 
with parallel walls during the passage of constant current through a nano channel 
junction. The basic idea in the derivation is depth averaging across the channel thick- 
ness, assuming that the electrical double layers are thin (e ^ 1) and confined to 
near-wall regions. Similar area-averaged models for nano channels with thin double 
layers have also been developed by l27l 67 82 . Recently, the thin-double-layer for- 



mulation for microchannels was (formally) extended to porous media by Mani and 
Bazant 49 by including the surface charge as a homogeneous background charge in 

In this limit, the (weakly) charged porous medium 
whose ion concentrations can be significantly de- 



the electroneutrality condition, 
behaves like a "bad membrane" 
pleted and enriched by the passage of current, since only a small fraction of the ions 
are involved in screening the surface charge. 
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A systematic extension of these results to general porous media requires more in- 



volved techniques to accurately treat the nonlinear terms in (3.6). At this point, the 
results from Section |3] are necessary to properly describe material transport in highly 
heterogeneous media, based on homogenization theory to derive effective Poisson- 



Nernst-Planck equations (3.6). A further advantage is that the system (3.6) is not 
restricted to a special geometry and rather can be adapted for general porous struc- 
tures by the tensors D. In the case of straight channels, the correction tensor D can 
be analytically obtained, as in Section [6j although for more complicated geometries, 
such as the irregular channels considered in Section [?[ the correction tensor D must 
be calculated numerically. 

In order to describe situations with thin electrical double layers compared to the 



mean pore size, we consider the limit e = / 1 ^ Q \n (3.6). In the general case of a 
polarizable solid matrix, one immediately sees that the limit e — > 0, does not reduce 
the complexity of the macroscopic formulation, when we apply this limit in the system 



( 5.1 ) and set e = in e(e, a). However, if we pass to the joint limit e, a — )■ 0, where the 
solid matrix is electrically insulating, then the porous media Poisson-Nernst-Planck 
system behaves almost like the classical PNP for e — 0, since we obtain the following 
leading order bulk approximation for salt density c, charge density p, and electric 
potential 0, i.e. 



= div (cDV(/>) , 



pdtc = div (dVc) - div ^^DV0j , (10-1) 

Q=pp + Ps- 



The first equation expresses charge conservation in the quasi-neutral bulk solution by 
setting the divergence of the current to zero. The second equation expresses total salt 
conservation. This description of bulk electrolytes with thin double layers is very well 
known and forms the basis for classical theories of electrochemical transport, based 
on the assumption of quasi-electroneutrality in the electrolyte p = 0, see Newman and 
Thomas- Alyea |57j. The third equation, however, is different and expresses quasi- 
electroneutrality of the entire porous composite, including not only the diffuse ionic 
charge p, but also the homogenized surface charge ps- Naively, one might expect such 
a relation to hold for thick double layers, but this is not necessarily the shown 
in the next section. 



Instead, the macroscopic electroneutrality condition. Equation (10.1), generally 



holds in the limit of thin double layers in charged porous media, as proposed by Mani 



and Bazant 49 . Physically, the reason is that the counter-ions screening the surface 
charge in a thin double layer provide an extra surface conductivity, proportional to 
the total diffuse double-layer charge, which is acted on by the same tangential electric 
field as in the nearby bulk solution. If the double layers were not thin, the electric 
field would be strongly perturbed by the diffuse charge throughout the pore, and the 
extra counter ions could not be viewed as simply providing extra conductivity for bulk 
ion transport. 

We briefly mention some simple cases of the thin-double-layer limit. In the case of 
straight channels from Section [6| the reduced porous media approximation becomes 
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in the two-dimensional case 

pdtC = pd'^^C- da:APsdxi(l)) , (10.2) 

Ps 

P = ■ 

P 

With the considerations in Section[6j a corresponding extension to the three-dimensional 
case is straightforward. Note that, although we have written the thin double-layer 



equations using the reference cell length and time scales 2.4 1, they have the same form 
at the macroscopic length and time scales [3^ 



11. Thick double layers at the pore scale . 

11.1. Membrane limit . Next we consider the limit of thick double layers, 
e = Ad/^ ^ 1, taken after the homogenization limit r — > 0. We first consider the 
"membrane limit" where the Debye length is much larger than the mean pore size, but 
much smaller than the macroscopic geometrical scale, e = er = \d/L <^ 1. Again, 
the upscaled equations take a simple form if the electric potential only exists in the 



solid phase, i.e., a 0. In that case, the macroscopic PNP equations (3.61 take the 
following form at leading order for thick double layers at the pore scale: 



(e^I)\7(j)^ ^p{c+ - c^) + Ps mQ. (11.1) 



pdfc^ = V- (dVc* ± c^I)V(j)) in n , (11.2) 



where we have rescaled to the macroscopic length and time scales (3.4 1. Asid e frorn 
the geometrical tensor D, this form of the macroscopic PNP equations 
reduces to a standard mathematical model for diffuse charge in membranes (15|17) and 
doped semiconductors |52 . As a result of the thick double layer limit over the pores, 
such materials act as "good membranes" and maintain roughly uniform conductivity 
outside the double layers, which are typically thin at the geometrical scale. 

11.2. Thin-film limit . For very thin porous films, we can also consider the 
limit where the double layer thickness is much larger than the macroscopic length 
scale, e = Xd/L ^ 1. For such thin films, Poisson's equation is replaced by 



(dV0 



(11.3) 

In this limit, the macroscopic electric field —Vcf), which drives electromigration of the 
ions, satisfies Poisson equation for an uncharged dielectric material with anisotropic 
permittivity. Mathematically, from Eq. ( |11.1[ ) we see that local electric field variations 

are very small, — 0{e ^), and show no significant effect from diffuse charge 

fiuctuations, which are corresponding small within the pores. It can be shown that a 



similar thick-double-layer, thin-film approximation (11.2) also holds true for the case 
where the electric potential is additionally defined in a polarizable solid matrix. 

A useful feature of this limit is that the potential satisfies Laplace's equation 



(11.2)2 (in transformed coordinates, as noted above), which is decoupled from the 
ion transport problem. The resulting electric field acts like an external potential fiow 
on the ions, leading to an advection-diffusion equation for the concentration profiles 
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( 11.2 )i. This is also a common approximation in gel electrophoresis, where it can also 



be justified for charged tracer particles outside thin macroscopic double layers (e <C 1), 
where the total conductivity is dominated by the fixed gel charge, according to the 



general thick double layer formulation (11.2). In two dimensions, this problem has 



an elegant mathematical structure in steady state, amenable to solution by conformal 
mapping [7 22 , which could be useful for this class of ion transport problems in 
charged porous media. 

12. Optimizing the conductivity for parallel straight channels. We are 

interested in the effective conductivity d'{x) of a binary symmetric electrolyte occu- 
pying the domain O with corresponding surface dfl. The domain Q is assumed to 
be a porous medium with porosity p. For simplicity, we first consider the pores to 
be straight (cylindrical) channels a nd t he case of an insulating porous matrix, that 
means D = M = e, see Section [s] |3.6| . For a current density J together with the 
electrostatic equations div J = and rot_E 0, where E — V0, we obtain 



div (.tV^) = , 



(12.1) 



the constitutive relation J — aE entered. The upscaled Poisson-Nernst-Planck equa- 
tions provide the current density J for a binary symmetric electrolyte, i.e. 



J := DVp + cDV0. 



(12.2) 



In order to determine in ( 12.1 ), we replace p in ( 12.2 1 by the Poisson equation 13.6^3 



with e^D instead of ikiie,a) as explained in Remark |2[ We obtain 
J = -DV ^div ^^DV0^ ^^^"^ ' 

P \ P J 



(12.3) 



The structure of equation ( 12.3 ) motivates to consider the eigenvalue problem for the 
Laplace operator, i.e. 



-AyUi{y) 
Ui{y) 



iUi{y) = 



in YP , 

on dYP n dV 



(12.4) 



We remark that it is not immediately clear what kind of boundary conditions are 
required in ([12.4 1. The boundary condition (12.4l2 has the advantage that it gives a 
lower bound |20|40| on the first eigenvalue 9i in ( 12.4 ) for the geometry defined by the 
pore phase Y^. We point out that instead of using the macroscopic Laplace operator 
Ax, we apply the microscopic Laplace operator Aj, — s'^A^ on the pore phase of the 
reference cell Y^. Hence, the eigenvalue 9 depends on the pore geometry which is 
the striking point for our optimization goal. Since the self-adjoint eigenvalue problem 



(12.4) is a regular Sturm-Liouville problem, we can use its solutions {u^}^ to generate 
an orthonormal basis in L^(57). Thus, for any function / e L'^{fl) we have 



(12.5) 



where equality is in the sense of L^. 
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Now, we can choose D as in Section [H] for straight channels, if we additionaUy 
assume that the electrostatic potential only exists in the electrolyte phase. Hence, 
after choosing / = 9x^0, the relation ( 12. 3[ ) becomes 



OO / 9 



(12.6) 



where equality holds again in the L -sense. We can now approximate (12.6) by only 



considering the first eigenvalue 6i of (12.4). That means, we obtain 

^2 



+ c {dx^(t),ui)ui ~ d,^^ps 



As a consequence, the conductivity cr can be approximated by 

(Til ~ P 



(12.7) 



(12.8) 



That means, optimization of the conductivity in direction of the straight pores is 
achieved by increasing e and di for given p, c and s. With the help of Cheeger's 
number h{flP), we obtain an additional tool to optimize the conductivity with respect 
to the geometry. In fact, it holds due to Cheeger 20 and Kawohl and Fridman 40 



> 



(12.9) 



Example 1: (Square) For a square Sa '■— [—a, a]^, Cheeger's number can be determined 



explicitely by h{Sa) = (4,2^^)0 • Moreover, we know that the first eigenvalue is 
^i('5'i) = 27r^. This indicates that the lower bound given by estimate (12.9) is not 
too sharp. However, it allows at least to obtain first insights for possible directions 

a, a] X [—6, 5], one immediately gets 



towards optimization of the conductivity (12.8) 
Example 2: (Rectangle) For a rectangle Ra. 



..b ■- 



the following Cheeger constant, see 41 



HRa.b) 



. + b- y/{a - 6)2 + nab 



(12.10) 



Hence, in order to optimize the conductivity (12.8) for a rectangle shaped pore Ra,b, 
we have to maximize h{Ra,b) what is equivalent to the minimization of a and b. If 
we assume that we are given a porous material of characteristic length b = I, then it 
immediately follows that h is maximal after minimizing the channel hight a > 0. 

13. Conclusion . We have applied a systematic, formal homogenization proce- 
dure for the Poisson-Nernst-Planck equations (3.1 ) for ion transport in charge porous 



media. The resulting upscaled macroscopic equations (3.6) have a similar form as the 
microscopic equations, except for two fundamental modifications: (i) The ionic diffu- 
sivities and mobilities, as well as the effective medium permittivity, become tensorial 
coefficients, which are explicitly connected to the microstructure by solving the peri- 
odic reference cell problem, and (ii) the total surface charge per volume appears as 
an extra "background charge" in the upscaled Poisson equation. The porous-medium 
PNP equations may find many applications in electrochemical and biological systems 
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involving ion transport in charged porous media, where effects of fluid flow can be 
neglected. Simplified equations for the limits of thin or thick double layers may also 
be appropriate in many cases. 

There are many interesting avenues for future work, building on these results. 
There is a substantial literature on rigorous bounds and approximations for the ef- 
fective diffusivity or conductivity of a composite medium 85^ , related to solutions of 
Laplace's equation with flux matching interfacial conditions. It would be challeng- 
ing and useful to derive analogous mathematical bounds and approximations for the 
effective diffusivities and mobilities of ions in a charged composite medium, which 
appear as tensorial coefficients in our porous-medium PNP equations. One might 
expect analogs of the Wiener bounds for anisotropic composites to hold for striped 
microstructures and analogs of the Hashin-Shtrikman bounds for isotropic microstruc- 
tures to hold for space-filling random sphere packings, although the appearance of an 
internal length scale for electrostatic interactions (the Debye screening length) com- 
plicates such simple geometrical constructions. 

It would also be valuable to find simple ways to approximate the solution to 
the reference-cell problem and thus derive simplified expressions for the tensorial 
diffusivities and mobilities. In the limit of thin double layers, this could be done using 
surface conservations laws, which are effective boundary conditions on the neutral 
solution obtained by singular perturbation methods 24 1. In the opposite limit of 
thick double layers, regular perturbation methods might be applied to capture effects 
of diffuse charge variations in the microstructure. 

We close by emphasizing the open challenge of deriving effective ion transport 
equations in more general situations using homogenization theory. We have already 
commented on the extension to concentrated solution theories based on the local den- 
sity approximation (for chemical interactions) and the mean-field approximation (for 
electrostatics) j9^. Going beyond these approximations in the microscopic equations 
can lead to non-local Nernst-Planck integral equations 33 34 or higher-order Poisson 
equations [TT|[2l|[35l[72 , whose upscaled form remains to be determined. Perhaps even 
more challenging, and more important for many applications, would be to predict the 
effects of fluid flow on the homogenized PNP equations, coupled to the Navier-Stokes 
equations with electrostatic body forces. 
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